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ON  THE  NYSTROM  DISCRETIZATION  OF  INTEGRAL  EQUATIONS  ON 
PLANAR  CURVES  WITH  CORNERS 

JAMES  BREMER 


Abstract.  The  Nystrom  method  can  produce  ill-conditioned  systems  of  linear  equations  and 
inaccurate  results  when  applied  to  integral  equations  on  domains  with  corners.  This  defect  can 
already  be  seen  in  the  simple  case  of  the  integral  equations  arising  from  the  Neumann  problem 
for  Laplace’s  equation.  We  explain  the  origin  of  this  instability  and  show  that  a  straightforward 
modification  to  the  Nystrom  scheme,  which  renders  it  mathematically  equivalent  to  Galerkin 
discretization,  corrects  the  difficulty  without  incurring  the  computational  penalty  associated  with 
Galerkin  methods.  We  also  present  the  results  of  numerical  experiments  showing  that  highly 
accurate  solutions  of  integral  equations  on  domains  with  corners  can  be  obtained,  irrespective 
of  whether  their  solutions  exhibit  bounded  or  unbounded  singularities,  assuming  that  proper 
discretizations  are  used. 


1.  Introduction 


It  is  well  known  that  Neumann  boundary  value  problems  for  Laplace’s  equation  can  be  formulated 
as  second  kind  integral  equations.  In  the  case  of  a  compact  simply-connected  planar  domain  f 1  with 
Lipschitz  boundary  dfl,  the  exterior  problem 


A u  =  0  in  flc 

du  on 

—  =  on  dQ 

ou 


(1.1) 


is  solvable  when  g  is  in  the  space  Lq(cAI)  of  square  integrable  functions  of  zero  mean  on  dil.  Here, 
v  is  the  outward-pointing  unit  normal  of  dfl  and  the  boundary  condition  must  be  understood  as  a 
limit  in  the  appropriate  non-tangential  sense.  If  the  additional  condition 


u(x )  =  O 


as  |x| 


(1.2) 


is  imposed,  then  the  solution  is  unique.  The  usual  integral  equation  method  exploits  the  observation 
that  the  solution  u  can  be  represented  uniquely  in  the  form  of  a  single-layer  charge  distribution  cr 
on  dil:  that  is,  as 


=  log  \x  —  y\  a(y)ds(y). 

^  Jd  n 


(1.3) 


The  charge  distribution  cr  which,  when  inserted  into  (1-3),  produces  the  solution  of  (1.1)  can  be 
determined  by  solving  the  integral  equation 


a(x)  +  Na(x)  =  2g(x),  x  £  Oft,  (1-4) 

where  N  is  the  operator 

Na(x)  =  -  (  (J^loglx -y\)  a{y)ds{y).  (1.5) 

nJdn\ot>x  J 

In  this  article,  we  will  consider  domains  dQ  which  are  piecewise  smooth  with  corner  points.  By 
corner  point,  we  mean  a  point  on  the  boundary  curve  dfl  at  which  the  limits  of  the  normal  derivative 
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as  one  approaches  from  the  clockwise  and  the  counter-clockwise  directions  exist  but  are  not  equal. 
For  these  domains,  the  kernel 

K{x,  y')  =  —  7C“—  log  \x  y\  (1.6) 

7T  UUx 

of  the  integral  operator  N  is  smooth  except  on  the  set  consisting  of  points  ( x ,  y)  £  dfl  x  df l  such 
that  either  x  or  y  is  a  corner  point. 

The  Nystrom  method  is  one  of  the  standard  approaches  to  the  discretization  of  the  integral 
equation  (1.4).  It  proceeds  by  fixing  an  appropriate  quadrature  formula  xi,...,  xn,  wi,...,  wn  and 
approximating  the  values  of  the  charge  distribution  a  at  the  quadrature  nodes  x±, . . . ,  xn  by  solving 
the  system 

n 

a(xi)  +y ^K(xj,Xj)a(xj)wj  =  2  g(xj),  i  =  (1.7) 

j= i 

of  n  linear  equations  in  the  n  unknowns  er( Xi), . . . ,  cr(xn). 

For  many  domains,  the  linear  system  (1.7)  is  an  effective  means  for  solving  the  integral  equa¬ 
tion  (1.4).  If  the  boundary  of  df l  is  twice  continuously  differentiable,  then  N  is  compact  as  an 
operator  on  the  Banach  space  C(dfl)  of  continuous  functions  on  dfl  endowed  with  the  uniform 
norm.  The  boundedness  and  invertibility  of  I  +  N  as  an  operator  on  C(dfl)  follow  from  the  Fred¬ 
holm  theory.  Moreover,  bounds  on  the  operator  norms  of  I  +  N  and  (J  +  N )_1  can  be  used  to 
estimate  the  Z°°  condition  number  of  the  linear  system  (1.7).  Convergence  results  follow  from  sim¬ 
ilar  considerations;  for  instance,  if  the  domain  fl  is  smooth,  then  exponential  convergence  can  be 
achieved  by  choosing  an  appropriate  quadrature  formula.  See,  for  instance,  [8]  for  a  discussion  of 
the  classical  Riesz  theory  and  its  application  to  the  numerical  solution  of  the  integral  equations  of 
potential  theory. 

When  the  boundary  dfl  is  merely  Lipschitz,  this  analysis  breaks  down.  In  this  case,  I  +  N  is 
bounded  and  invertible  as  an  operator 

I  +  N  :  L20(dfl)  ^  L20(dfl)-, 

see,  for  instance,  [5].  But  now  (/  +  TV)-1  is  not  bounded  with  respect  to  the  uniform  norm; 
indeed,  there  are  functions  g  £  C(dfl)  such  that  the  charge  distribution  a  satisfying  the  integral 
equation  (1.4)  is  unbounded  (see  Section  4.3).  So  the  usual  estimates  on  the  l°°  condition  number 
of  the  linear  system  (1.7)  no  longer  apply.  This  difficulty  led  to  conditioning  problems  in  [3]  and 
it  appears  to  be  the  motivation  for  [4],  which  notes  that  the  Nystrom  method,  when  applied  to 
boundary  integral  operators  on  domains  with  corners,  can  produce  ill-conditioned  linear  systems 
and  inaccurate  results. 

Stability  can  be  restored  by  discretizing  I  +  N  as  an  operator  on  L2(dfl)  rather  than  C(dfl).  This 
can  be  accomplished  by  switching  from  Nystrom  to  Galerkin  discretization.  That  is,  by  choosing  an 
appropriate  orthonormal  basis  il>i ,  •  •  • , tpn  and  discretizing  the  integral  operator  I  +  N  via  the  nxn 
matrix  A  with  entries 

Aij  =  Sij+  /  K(x,  y)ipj(y)il>i(x)ds(y)ds(x) . 

J  d£l  J  dfl 

Bounds  for  the  condition  number  of  A  as  a  function  of  singular  values  of  the  operator 

I  +  N  :  L2{dfl)  ->  L2(dfl) 

can  be  readily  obtained.  There  is,  however,  a  disadvantage  to  this  approach:  Galerkin  discretizations 
are  computationally  more  expensive  and  more  complicated  to  construct  than  Nystrom  discretizations 
because  of  the  need  to  evaluate  double  integrals  numerically. 

It  is  the  purpose  of  this  article  to  show  that  a  second  remedy,  which  does  not  suffer  from  this 
disadvantage,  is  available.  In  particular,  the  Nystrom  scheme  can  be  modified  to  produce  discretiza¬ 
tions  which  reflect  the  properties  of  I  +  N  as  an  operator  on  Lg(<9£4).  The  modified  scheme  results 
in  well-conditioned  linear  systems  approximating  (1.4)  on  domains  with  corners  whose  solutions  are 
highly  accurate.  Moreover,  this  method  is  simpler  and  more  efficient  than  Galerkin  discretization 
in  that  it  does  not  require  the  numerical  evaluation  of  double  integrals. 
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There  are  further  difficulties  complicating  the  discretization  of  integral  equations  on  domains 
with  corners:  the  operator  N  is  not  compact  and  both  the  kernel  K{x,y)  and  solutions  of  the 
equation  (1.4)  are  singular.  It  has  long  been  standard  practice  to  address  the  first  problem  by 
forming  a  modified  boundary  curve  dfls  which  excludes  from  dfl  a  (5-neighborhoocl  of  each  corner 
point,  introducing  the  compact  integral  operator 

Ns(x)  =  -[  f^-log\x-y\]a(y)ds(y),  (1.8) 

7T  JdQs  \OVx  > 

and  considering,  in  lieu  of  (1-4),  the  integral  equation 

I  +  Ns  =  2 g(x),  x  €  dfls-  (1.9) 

It  is  well  established  that  solutions  as  of  (1.9)  converge  to  those  of  (1.4)  in  a  distributional  sense; 
that  is, 

/  log  | x  -y | as (y) ds (y)  ->  /  log  \x  -  y\a(y)ds(y)  as  6  ->  0 
J dQs  J 9Q 

for  x  £  flc.  To  combat  the  difficulties  associated  with  singular  solutions  and  kernels,  the  Nystrom 
discretization  of  (1.9)  is  usually  formed  using  graded  mesh  quadratures;  that  is,  with  the  help  of 
quadratures  of  various  types  which  become  increasingly  dense  near  corner  points. 

Moreover,  while  the  discretization  of  boundary  integral  equations  on  domains  with  corners  via 
the  Nystrom  method  and  graded  mesh  quadrature  can  be  made  accurate  and  stable,  this  is  not  by 
itself  a  reasonable  approach  to  the  solution  of  boundary  integral  equations  on  large-scale  domains 
with  corners.  The  systems  of  linear  equations  arising  from  such  domains  are  excessively  large  and 
must  be  compressed  in  some  fashion.  In  Section  3,  we  present  a  scheme  for  the  compression  of 
linear  systems  obtained  by  discretizing  boundary  integral  equations  on  domains  with  corners  which 
is  based  on  ideas  introduced  in  [2]  and  [3].  In  conjunction  with  fast  multipole  methods  and  fast 
direct  solvers,  it  allows  for  the  very  rapid  solution  of  boundary  integral  equations  on  large-scale 
domains  with  corners.  Although  we  only  consider  the  Neumann  problem  for  Laplace’s  equation 
here,  the  scheme  of  Section  3  applies  to  a  wide  class  of  boundary  value  problems;  see,  for  instance, 
[3]  where  a  similar  approach  is  used  to  solve  the  Dirichlet  problem  for  the  Helmholtz  equation  at 
low  wavenumbers.  Moreover,  certain  components  of  these  schemes  generalize  readily  to  the  three- 
dimensional  environment.  The  specific  approach  discussed  in  Section  3  can  be  readily  extended  to 
this  setting,  for  instance. 

This  article  is  structured  as  follows.  In  Section  2,  we  discuss  the  Nystrom  discretization  of  integral 
operators  on  domains  with  corners.  Section  3  gives  the  details  of  a  scheme  for  the  compression  of 
systems  of  linear  equations  arising  from  boundary  integral  equations  on  domains  with  corners.  In 
Section  4,  we  present  the  results  of  numerical  experiments  showing  that  the  boundary  integral 
equation  (1.4)  can  be  solved  stably  and  to  high  accuracy  with  graded  mesh  quadrature  provided 
appropriate  discretizations  are  used.  Finally,  we  close  with  a  few  remarks  on  the  numerical  solution 
of  integral  equations  on  domains  with  corners  in  Section  5. 

2.  Discretization  of  integral  equations  on  domains  with  corners 

2.1.  Discretization  of  compact  operators  acting  on  L2  spaces.  If  (X,y)  is  a  measure  space, 
then  we  say  that  a  mapping  4>  of  M  C  L2  ( X ,  dy)  into  C”  preserves  inner  products  provided 

$(/)  '  ®(9)  =  [  f{x)g{x)dy{x) 

Jx 

for  all  f,g£  M. 

Consider  a  compact  operator  S  :  L2(X,  dy)  — >  L2(X,  dy).  It  is  well  known  that  S  admits  singular 
value  decomposition.  We  will  denote  by  {Xj}  the  singular  values  of  S  and  for  each  j,  we  will  let  Uj 
and  Vj  be  the  left  and  right  singular  functions  corresponding  to  Xj.  In  the  event  that  the  rank  of  S 
is  finite,  we  set  Xj  =  0  and  Uj  =  Vj  =  0  for  all  j  greater  than  the  rank  of  S.  Moreover,  we  will  let 
Ue  and  Ve  be  the  subspaces  spanned  by  {uj  \  Xj  <  e}  and  [vj  \  Xj  <  e},  respectively. 
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We  say  a  matrix  A  in  Cnx™  is  an  inner  product  preserving  discretization  of  S  of  accuracy  e  if 
there  exist  inner  product  preserving  mappings  $  of  U£  into  Cn  and  T  of  V£  into  Cn  such  that  the 
diagram 

Ve  C  L2(X,g)  U£  C  L2{X ,n) 

Cn  - C" 

commutes.  This  is  equivalent  to  requiring  that 

AA>{f)  •  $(5)  =  [  Sf(x)g{x)dn(x) 

Jx 

for  all  /  £  V£  and  g  £  Ue.  One  of  the  key  properties  of  such  a  discretization  is  that 

A*AV(f)  ■  T(.g)  =  A*(f)  ■  A*(g)  =  f  Sf(x)Sg(x)dfx(x )  =  [  S*Sf(x)g(x)dg(x ) 

■lx  Jx 

for  all  f,g  £  V£,  from  which  it  follows  that  the  singular  values  of  the  restriction  of  the  matrix  A  to 
'li(Ve)  are  those  of  the  restriction  of  the  operator  S  to  V£. 

Remark  2.1.  Galerkin  discretization  is  perhaps  the  easiest  way  to  produce  an  inner  product  pre¬ 
serving  discretization  of  a  compact  operator  S  :  L2{X,dg)  — >  L2{X,dy).  Let  /i,...,/n  be  an 
orthonormal  basis  spanning  Ve  and  let  gi , . . . ,  gn  be  an  orthonormal  basis  spanning  Ue .  The  map¬ 
ping  which  takes  functions  in  Ve  to  their  coefficient  expansions  with  respect  to  the  basis  {fj}  is  inner 
product  preserving,  as  is  the  mapping  which  takes  functions  in  Ue  to  their  coefficient  expansions  in 
the  basis  {gi}.  The  Galerkin  discretization  of  S  with  respect  to  these  bases,  which  is  the  nxn  matrix 
A  with  entries 

Aij  =  /  Sfi{x)gj{x)dp{x), 

Jx 

is  then  an  inner  product  preserving  discretization  of  S  of  accuracy  e. 


$ 


2.2.  Inner  product  preserving  Nystrom  discretizations.  Now  suppose  that  S'  is  a  compact 
integral  operator  L2{X,dp.)  — >  L2[X,dfi)  whose  kernel  L(x,y)  is  defined  pointwise  everywhere  and 
let  Ve  and  Ue  be  as  in  the  preceding  section.  A  quadrature  rule  of  the  form 


/  f{y)dn{y)  « 

lx 


(2.1) 


with  positive  weights  induces  a  mapping  $  of  L2{X,dg)  into  Cn;  to  wit, 


*(/) 


(  f{x  I)y/Wf  \ 

f(X2)V^ 


\  f{xn)y/wh  ) 


If  the  quadrature  rule  (2.1)  is  exact  for  products  of  the  form  fg  with  f,g£Ve  and  hf  with  h,  f  £  U£, 
then  <I>  is  an  inner  product  preserving  mapping  of  U£  and  V£  into  C". 

Moreover,  the  choice  of  quadrature  Xi, . . . ,  xn,  Wi, . . . ,  wn  induces  a  mapping  of  the  space  of 
integrals  operators  on  L2{X,dg)  with  pointwise  defined  kernels  into  the  space  of  operators  on  Cn. 
In  particular,  S  is  mapped  to  the  nxn  matrix  A  with  entries 


Aij  L{i} 

If  the  quadrature  rule  is  sufficiently  dense  to  integrate  the  products  of  the  kernel  with  the  functions 
in  Ue  and  V£,  then  we  will  have 


A®  (f)  ■  $g{x)  =  [  Sf(x)g(x)ds(x) 

Jx 

for  f  £  V£  and  g  £  XJ£,  that  is,  A  will  be  an  inner  product  preserving  discretization  of  the  integral 
operator  S  of  accuracy  e. 
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This  notion  of  discretization,  when  applied  to  the  operator  I  +  N,  leads  to  the  n  x  n  matrix  B 
with  entries 

Bij  —  T  K (/X'l ,  Xj  )  yj Wj ,  (2.2) 

where  K(x,  y)  is  the  kernel  (1.6)  of  N.  This  is  in  contrast  to  the  standard  Nystrom  discretization 
of  I  +  N,  which  is  the  matrix  C  with  entries 

Cij  —  Sij  T  K(xi,Xj)wj.  (2.3) 

The  salient  difference  between  the  matrices  B  and  C  is  that  the  singular  values  of  B  approximate 
those  of  I  +  N  while  the  singular  values  of  C  do  not.  Since  the  operator  I  +  N  and  its  inverse 
have  well-behaved  L2(dfl)  operator  norms  on  most  curves  of  interest,  it  follows  that  the  matrix  B  is 
well-conditioned  in  typical  cases,  including  on  many  singular  domains.  The  matrix  C,  on  the  other 
hand,  is  generally  ill-conditioned  when  dLl  is  singular. 

The  solution  of  the  integral  equation  (1.4)  using  the  inner  product  preserving  discretization  B 
proceeds  by  letting  y  be  the  vector  of  length  n  with  entries 

Vi  =  g{xi)y/wi 

and  inverting  the  linear  system  Bz  =  y.  The  entries  of  the  resulting  vector  z  are  the  values  of  the 
charge  distribution  er  satisfying  (1.4)  at  the  quadrature  nodes  Xi  scaled  by  the  square  roots  of  the 
corresponding  quadrature  weights;  that  is, 


Zi  =  a(xi)y/wi~. 

In  other  words,  the  analog  of  (1.7)  is  the  system 

n 

<r{xj)y/wl  +  ^K(xi,xj)y/wfy/w]a(xj)  =  2  gjx^y/wf,  i  =  l,...,n,  (2.4) 

j= i 

of  n  linear  equations  in  the  n  unknowns 

<T{xi)y/wi,  Cr(x2)y/W2,  ••.,  <7{xn)^/Wn. 

We  will  refer  to  linear  systems  of  this  type  as  inner  product  preserving  discretizations  of  the  integral 
equation  (1.4).  Note  that  the  vector  obtained  from  solving  the  linear  system  (2.4)  is  the  image  of 
the  charge  distribution  cr  under  the  inner  product  preserving  mapping  4>;  i.e.,  we  solve  for  a  in  an 
L2(dn)  sense  rather  than  in  a  pointwise  sense. 

Remark  2.2.  We  have  neglected  certain  details  in  this  section  in  order  to  simply  the  discussion. 
Specifically,  in  finite  precision  arithmetic  it  is  inevitable  that  the  quadrature  formula  (2.1)  will  not 
hold  exactly  for  the  requisite  integrands.  This  means  that  our  discretizations  will  suffer  from  inac¬ 
curacies  both  as  a  result  of  finite  rank  approximation  and  from  quadrature  error.  Limitations  of  this 
type,  however,  apply  to  all  methods  for  the  discretization  of  integral  operators. 

Remark  2.3.  Inner  product  preserving  Nystrom  discretizations  can  be  produced  in  the  case  of 
integral  operators  whose  kernels  are  not  defined  pointwise  provided  the  appropriate  interpolatory 
quadratures  are  available. 


3.  Compression  of  linear  systems  arising  from  integral  equations  on  domains  with 

CORNERS 

We  now  describe  a  mechanism  for  the  compression  of  systems  of  linear  equations  arising  from 
boundary  integral  equations  on  domains  with  corners.  It  is  a  variant  of  the  schemes  of  [3]  and  [2] 
which  operates  via  two  tools,  generalized  quadrature  and  charge  bases.  Generalized  quadrature  and 
charge  bases  are  discussed  in  Sections  3.1  and  3.2,  respectively,  and  the  compression  scheme  proper 
is  detailed  in  Section  3.3. 
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3.1.  Generalized  quadrature.  We  define  the  e-rank  of  an  nx  m  matrix  A  to  be  the  least  integer 
k  such  that  the  k  largest  singular  values  of  A  are  greater  than  e,  assuming  that  such  an  integer 
exists.  If  it  does  not,  then  we  say  that  the  e-rank  of  A  is  min(m,n). 

If  /i , . . . ,  fm  is  a  collection  of  functions  in  L2  (X,  p)  and  <I>  is  an  inner  product  preserving  mapping 
of  the  span  of  the  fj  into  C",  then  we  define  the  e-rank  of  /i, . . . ,  fm  to  be  the  e-rank  of  the  n  x  m 
matrix 

A  =  (  Hfl)  ■  ■  ■  Hfm)  ) 

whose  columns  consist  of  the  images  of  the  fj  under  the  embedding  $. 

As  described  in  the  last  section,  a  quadrature  formula  Xi, . . . ,  xn,  w\, . . . ,  wn  with  positive  weights 
which  integrates  products  of  the  functions  fi , . . . ,  fm  induces  an  inner  product  preserving  mapping 
of  the  span  of  the  fj  into  C" .  The  following  theorem  on  the  existence  of  efficient  quadrature  formulas 
is  an  almost  immediate  consequence  of  that  fact. 


THEOREM  3.1.  If  x  1, . . . ,  xn,  wi, . . . ,  wn  is  a  quadrature  rule  with  positive  weights  which  inte¬ 
grates  products  of  the  functions  /i ,  •  •  • ,  fm  G  L2  ( A,  dpi) ,  then  for  every  e  >  0  there  exists  a  quadrature 
rule  y  1, . . . ,  yk,  V\, . . . ,  Vk,  whose  length  k  is  the  e-rank  of  the  collection  /i, . . . ,  fm,  such  that 

i  1 2 


E 


lx 


f{x)dy{x)  -  f(Vj)vj 
j= i 


<  (1  +  mk{m.  —  k))  e2 . 


Quadratures  of  this  type  can  be  constructed  by  finding  a  sparse  least  squares  solution  to  a  system 
of  linear  equations.  Let 


f2{xi)^/wi 


fl{x2)^/W2 

f2{x2)^W2 


fl(xn)^Wn  \ 
f2{xn)\/Wn 


\  fm{Xl)^/wf  fm{x2)^/wf 


,/  m  (x  n  )  \/V  'n  J 


(3.1) 


and 

bi 
b2 

\bm  ) 

where 

bj  =  /  fj(x)dfi(x),  j  =  1, . . . ,  to. 

J  x 

If  the  e-rank  of  the  collection  f\ , . . . ,  fm  is  k,  then  a  vector  2  G  Cn  with  no  more  than  k  nonzero 
entries  such  that 

\\Fz  —  b\\2  <  e\J\  +  mk(m  —  k) 

can  be  constructed  stably  in  0(nm2)  floating-point  operations  (see  [6]).  If  we  let  i i, . . . ,  A  denote 
the  indices  of  the  nonzero  entries  of  the  vector  2  and  set 


Vj  =  xij  ,  j  =  l,...,fc, 

and 

Vj  =  zij^w~,  j  =  l,...,k, 

then  j/i, . . . ,  j/fc,  Ui, . . . ,  Ufe  is  a  fc-point  quadrature  for  the  functions  /i, . . . ,  fm.  See  [1]  for  an  extended 
discussion  of  the  construction  of  such  quadratures. 

Remark  3.1.  Note  the  scaling  by  square  roots  of  quadrature  weights  in  (3.1).  As  in  the  case 
of  integral  operators,  the  use  of  inner  product  preserving  embeddings  of  this  type  leads  to  bounds 
on  the  condition  number  of  the  matrix  (3.1).  Their  omission,  by  contrast,  generally  results  in  ill- 
conditioning.  For  instance,  if  the  functions  fj  are  taken  to  be  the  monomials  x on  the  interval 
[—1,1]  and  square  roots  are  omitted,  then  (3.1)  becomes  a  Vandermonde  matrix. 
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3.2.  Charge  bases.  In  this  section,  we  show  that  a  basis  spanning  the  space  of  restrictions  of 
solutions  of  the  integral  equation  (1.4)  to  a  corner  region  can  be  constructed  under  a  mild  assumption 
on  the  right-hand  sides  of  (1.4).  We  will  refer  to  a  basis  of  this  type  as  a  charge  basis  for  the  corner 
region. 

Let  To  be  a  corner  region  which  is  part  of  the  boundary  curve  9f l  and  let  B  be  the  disc  of 
minimum  radius  centered  at  the  corner  point  of  To  which  contains  To-  Denote  by  Id  the  portion  of 
the  contour  9fi  contained  in  2 B  \  B ,  and  by  T2  the  portion  of  90  contained  in  the  complement  of 
the  disc  2 B.  Figure  1  depicts  the  situation.  Our  assumption  on  right-hand  side  g(x)  of  (1.4)  is  that 
it  is  the  normal  derivative  of  a  function  satisfying  the  Laplace  equation  in  the  disc  2 B. 


Figure  1:  The  geometry  assumed  for  the  charge  basis  construction. 


The  boundary  integral  equation 

a(x)+  f  K(x,y)a(y)ds(y)  =  2g(x),  x  €  90,  (3.2) 

JdQ. 

can  be  rearranged  as 

cr{x)+  f  K(x,y)a(y)ds(y)  =  2g(x)~  f  K(x,y)a(y)ds(y)~  f  K(x,y)a(y)ds(y),  x  €  90.  (3.3) 

1^2  -/Ti 

The  charge  basis  construction  procedure  operates  by  considering  the  restricted  integral  equation 


cr(x)+  K(x,y)a(y)ds(y)  =  2g(x)-  K(x,y)a(y)ds(y)~  K(x,y)a(y)ds(y),  ier0,  (3.4) 
jTq  ^  r  2  JTx 

and  introducing  series  approximations  for  the  terms  on  the  right-hand  side.  By  solving  this  restricted 
equation  for  each  term  of  each  approximation,  we  form  a  set  of  functions  which  spans  the  space  of 
solutions  of  equation  (3.4).  Of  course,  the  restriction  of  any  solution  of  the  integral  equation  (3.2) 
to  r0  satisfies  (3.4),  so  such  a  collection  of  functions  is  a  charge  basis. 

The  first  approximation  depends  on  our  assumption  on  the  right-hand  side,  namely,  that  g  is  the 
restriction  to  dfl  of  the  normal  derivative  of  a  function  u  which  is  harmonic  in  the  disc  2 B.  It  is 
well  known  that  u  can  be  represented  as  a  multipole  expansion 

OO 

ti(x)  =  ^  oyr-'  cos (jO)  +  PjV-’  sin (jO) 

J=o 

in  the  disc  2 B.  Here,  (r,  6)  denotes  the  usual  polar  coordinate  system  with  respect  to  the  center  of 
the  disc  B.  For  points  x  in  B,  this  expansion  achieves  exponential  convergence;  in  particular,  there 
is  a  constant  C  such  that  for  x  £  B  we  have 


N 

u(x )  —  ayr7  cos  (jO)  +  PjT-’  sin  (jO) 

l=o 


<  C2~N . 
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It  follows  that  g  can  be  approximated  on  To  as  the  normal  derivative  of  a  finite  sum  of  multipoles; 
that  is,  as 

Nl  8  8 
E  a0  E  C0SU8))  +Pj  (V  sin (jO))  , 

j= i 


where  the  number  of  terms  Ni  grows  logarithmically  with  the  precision  of  the  approximation. 

Similar  considerations  apply  to  the  second  term  in  equation  (3.4).  Since  K(x,y)  is  harmonic 
and  To  is  separated  by  a  distance  equal  to  its  radius  from  F 2,  for  x  in  To  we  can  introduce  the 
approximation 


>r2 


AT, 

K(x,y)a(y)dy  «  E  (y^r7  cos  (jO)  +  yjr^  sin(j0)) 
1=0 


(3.5) 


where  once  again  (r,6)  is  the  polar  coordinate  system  with  respect  to  the  center  of  the  disc  B.  The 
number  of  terms  _/V2  in  (3.5)  also  grows  as  O  (log2  (e))  with  the  desired  precision  e. 

The  approximation  of  the  last  term  in  (3.4)  requires  a  different  strategy  since  T0  is  not  separated 
from  Ti.  Here  we  exploit  the  fact  that  the  integral  operator  L2{Ti)  — »  L2(Fq)  defined  by 


Rcr(x)  =  f  K(x,y)a(y)ds(y)  (3.6) 

Jr1 

is  compact  and,  barring  complex  geometry,  of  low  rank.  This  means  that  the  operator  R  can  be 
approximated  via  a  quadrature  formula 


K(x,y)a(y)ds(y) 


n3 

EA^  1  yj)a{yj)wj- 
i= 1 


(3.7) 


The  number  of  terms  in  (3.7)  depends  on  the  desired  accuracy  of  the  approximation  and  the  decay 
of  the  singular  values  of  integral  operator  R.  In  most  cases  of  interest,  N3  grows  as  O  (logg(e)), 
where  e  is  the  accuracy  achieved;  see  [3]  for  a  detailed  estimate. 

It  follows  that  a  charge  basis  can  be  formed  by  solving  the  restricted  integral  equation  (3.4)  for 


the  normal  derivatives 

d 

— rJ  cos  (id) 
dis  v  ’ 

d  ■ 

and  —  r3  sin(j6»),  j  =  1. . .  .V,, 

(3.8) 

the  multipoles 

r7  cos  (jO) 

and  r7  sin(j'0),  j  =  0, . . . ,  N2, 

(3.9) 

and  the  kernel  functions 

K(x,yj),  yi,...,yN3 

(3.10) 

It  is  easy  to  see  from  our  earlier  estimates  that,  in  typical  cases,  the  dimension  of  the  resulting  basis 
is  O  (logg  (e)),  where  e  is  the  desired  precision  for  the  basis. 


3.3.  The  compression  procedure.  The  mechanisms  discussed  in  the  two  preceding  sections, 
charge  bases  and  generalized  quadrature,  can  be  used  to  compress  the  systems  of  linear  equa¬ 
tions  resulting  from  boundary  integral  equations  on  domains  with  corners.  We  now  describe  one 
algorithm  for  doing  so;  it  is  a  variant  of  that  introduced  in  [2]  and  amounts  to  a  local  direct  solver 
for  corner  regions. 

Setup  for  the  procedure 

We  will  consider  the  discretization  of  the  integral  equation  (1.4)  on  a  contour  dfl  which  is  piecewise 
smooth  with  a  corner  region  To-  By  corner  region  we  mean  a  compact  connected  subset  of  dfl  which 
contains  a  corner  point  in  its  interior.  Let  r  :  [—5,  <5]  — >  Tq  be  a  parameterization  of  the  corner  region 
r0.  Also,  let  X\, . . . ,  xn,  Wi, . . . ,  wn  be  a  quadrature  rule  for  functions  on  the  corner  region  r0  and 
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let  xn+i, . . . ,  Xn+m ,  wn+ 1, . . . ,  Wn+rn  be  a  quadrature  rule  for  functions  on  the  remaining  portion  of 
the  curve  d SI  \  r0,  so  that  we  have  the  approximation  formulae 


and 


f{x)ds{x)  «  f(xj)wj 
i=i 


/9n\r0 


f(x)ds(x ) 


n+m 

j—n+l 


These  quadratures  lead  to  the  system 


n+m 

cr(xj)^/wl+  Y  K(xi,yj)a(xj)^/wis/wj  =  2 g(xj)^/wi,  i=l,...,n  +  m,  (3.11) 

j— i 

of  n  +  m  linear  equations  in  n  +  m  unknowns,  which  discretizes  the  integral  equation  (1.4). 

The  purpose  of  this  algorithm  is  to  replace  the  linear  system  (3.11)  with  a  compressed  system 
of  l  +  m  linear  equations  in  l  +  in  variables,  where  l  <C  n.  This  is  accomplished  by  constructing  a 
quadrature  rule  t/i, . . . ,  yi,  Vi, . . . ,  vt  such  that  the  restrictions  of  solutions  a  of  the  equation  (1.4)  to 
T o  are  fully  determined  by  their  values  at  the  quadrature  nodes  yi,...,yi-  This  quadrature  is  used 
to  discretize  the  self-interaction  of  T0  as  an  l  x  l  matrix  and  the  interactions  of  T0  with  dSl  \  T0 
through  l  x  to  and  to  x  l  matrices. 

Step  one:  construction  of  a  charge  basis 

The  restricted  integral  equation  (3.4)  on  To  is  discretized  using  the  quadrature  rule  with  nodes 
x\, . . . ,  xn  and  weights  W\, . . . ,  wn.  The  resulting  linear  system  is  then  solved  for  the  functions  (3.8), 
(3.9),  and  (3.10).  The  parameters  N\,  N-2,  and  JV3,  which  control  the  number  of  right-hand  sides, 
can  be  taken  to  be  quite  small  because  of  the  exponential  convergence  of  the  various  approximations 
of  the  preceding  section.  For  the  experiments  of  this  article,  Ni,  N2,  and  N2j  were  all  taken  to  be 
10.  As  discussed  in  Section  3.2,  the  resulting  solutions,  which  we  will  denote  by  <j\,  cr2, . . . ,  crk,  form 
a  charge  basis  for  the  integral  equation  (1.4)  over  the  region  To- 

Step  two:  formation  of  a  generalized  quadrature 

In  this  step,  a  quadrature  rule  for  integrals  of  the  form 


'To 


K(x,  y)cr j(y)ds(y),  j  =  l,...,k, 


x  G  dfi  \  r0, 


is  constructed.  To  accomplish  this,  we  take  advantage  of  the  fact  that  the  function  /  :  To  — ►  K. 
defined  by  f(y)  =  K{x,y)  is  piecewise  smooth  with  a  single  discontinuity  at  the  corner  point  of  To 
when  x  £  dfl  \  To-  This  means  that  f(y)  can  be  approximated  on  To  by  the  images  of  piecewise 
polynomials  on  the  intervals  [ — A,  0]  and  [0,<S]  under  the  parameterization  r  :  [— (5,  <5]  — ►  To-  In  the 
first  instance,  an  orthonormal  basis  (f>\ ,...+;  for  the  space  spanned  by  the  functions 

aj(x)p(x)  +  q(x),  j  =  l,..,,k,  (3.12) 

where  p{x)  and  q{x)  are  images  of  piecewise  polynomials  of  a  fixed  order,  is  formed.  Then,  the 
procedure  described  in  Section  3.1  is  applied  to  the  functions  in  order  to  generate  the 

desired  quadrature  formula  yi, . . . ,  yi,  v\, . . . ,  vp  This  quadrature  rule  is  generally  much  smaller 
than  the  initial  quadrature  Xi, . . . ,  xn,  Wi, . . . ,  wn.  For  the  calculations  of  this  article,  we  took  the 
degree  of  the  polynomials  to  be  9. 

Step  three:  construction  of  two  transformation  matrices 
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Let  $  denote  the  inner  product  preserving  mapping  of  the  span  of  the  (f>i , . . . ,  (f>i  into  Cl  defined 


by 


(  ai\ 

0.2 


$  = 

V  «I  ) 

and  let  T  be  the  mapping  of  the  span  of  0i , . . . ,  (j>i  into  Cl  induced  by  the  quadrature  y±, . . . ,  yi,  v\, . . . , 


*(/)  = 


(  /(yi)\Ah  \ 

f{y2)y/V2 


V  f(yi)V^i  / 


The  mapping  T  is  not  inner  product  preserving.  Nonetheless,  as  described  in  [3],  the  matrices 
and  <f>'f'_1  exist  and  are  well-conditioned.  The  matrix  Td)'1  :  Cl  — >  Cl  takes  the  coefficient 
expansion  of  a  function  /  with  respect  to  the  basis  <j>i , . . . ,  fa  to  its  scaled  values  at  the  quadrature 
nodes  (here,  by  scaled  values,  we  mean  the  values  of  the  function  /  at  the  quadrature  nodes  yi, . . .  ,yi 
scaled  by  the  square  roots  of  the  corresponding  quadrature  weights).  Similarly,  takes  the 

scaled  values  of  a  function  /  in  the  span  of  the  (j>i ,  ■  •  • ,  <t>i  to  its  coefficient  expansion  with  respect 
to  01,.  ,  (j>i. 

In  this  step,  we  construct  the  l  x  l  matrices  and  by  first  forming  the  matrix 


T'h"1 


/  Myi)V^  <h(yi)V uT 

Mv2)V^ 


Myi)Vn  \ 

MV  2)y/V2 
Mvi)\Mi  ) 


V  0i (yi)y/vi  02 {yi)Mh  • 

and  then  inverting  it  to  obtain  ‘hfE'-1. 

Step  four:  formation  of  the  compressed  linear  system 

In  this  final  step,  the  system  (3.11)  of  n  +  m  linear  equations  in  n  +  m  unknowns  is  replaced  by 
a  compressed  system  of  l  +  in  linear  equations  in  l  +  m  unknowns  of  the  form 


si 

S2 


Mi  M2 

Mi  M2 


51 

52 


=  2 


ti 

t'2 


(3.13) 


where: 

Aip  is  an  l  x  l  matrix  discretizing  the  operator  Tn  :  L2(To)  — >  L2(To)  defined  by 

Tucr(x)  =  f  K(x,y)a(y)ds(y)\ 

•'T0 

A 12  is  an  l  x  m  matrix  discretizing  the  operator  T12  :  L2(dQ  \  To)  — >  L2(To)  defined  by 

Ti2<j(x)  =  [  K(x,y)a(y)ds(y ); 

Jo  n\r0 

A‘2[  is  an  to  x  l  matrix  discretizing  the  operator  T2i  :  L2(T0)  — »■  L2(90  \  T0)  defined  by 

T21a(x)=  /  K(x,y)a(y)ds(y);  and 
Jr0 

A22  is  an  m  x  to  matrix  discretizing  the  operator  T22  :  L2(D \  To)  — >  L2(dfl  \  To)  defined  by 

T22<j(x)  =  [  K(x,y)a(y)ds(y). 

Jan\r0 
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The  l  x  l  matrix  An  is  constructed  by  forming  the  product 

An  =  B  ■(<$>■$ -1)  ,  (3.14) 

where  is  the  transformation  matrix  from  the  preceding  step  and  B  is  the  l  x  l  matrix  whose 

entries  are 

Bij  =  sfvi  /  K(y.i,  x)cj>j(x)ds(x),  i  =  l,...,l,  j  =  1, . . . ,  l. 

Jr0 

Note  that  the  entries  of  B  can  be  calculated  with  the  initial  quadrature  rule  x\ , . . .  ,xn,wi, . . . ,  wn 
for  the  corner  region  IV  That  An  approximates  the  operator  Tn  can  be  seen  easily.  The  matrix 
$'I'~1  maps  the  scaled  values  of  functions  in  the  span  of  the  basis  <j>\, . . .  ,<f>i  to  their  coefficient 
expansions  and  the  matrix  B  takes  coefficient  expansions  to  the  scaled  values  of  their  images  under 
the  operator  Tn- 

The  matrix  A22,  which  corresponds  to  the  self- interaction  of  the  curve  segment  <911  \  To,  is  formed 
in  the  usual  fashion;  that  is,  A 22  is  the  m  x  m  matrix  with  entries 

(^22  )ij  =  K(xn+i,xn+j)y/wn+iy/wn+j,  i  =  l,...,m,  j  =  l,...m. 

The  matrices  An  and  A-2n  which  approximate  the  interactions  of  To  with  (9n\To,  are  discretized 
through  the  quadrature  yi, . . .  yi,  Vi, . . . ,  Vi-  Specifically,  An  is  taken  to  be  the  l  x  m  matrix  with 
entries 

(An)i:j  =  K(yi,xn+j)y/wn+js/vl,  i  =  j  = 

and  A21  is  taken  to  be  the  m  x  l  matrix  with  entries 

(^21)^-  =  K{xn+i,yj)^/wn+i^/v~j,  i  —  1 , . . . ,  to,  j  =  l,...,l. 

Finally,  we  form  the  right-hand  side  of  the  system  (3.13)  by  setting 


(  g(yi)\Bh  \ 

/  g(xn+i)y/wn+1  \ 

t\  — 

5(2/2)^ 

11 

£ 

c5 

g{xn+2)s/Wn+2 

v  g(yi)V vi  ) 

The  unknowns  in  the  compressed  linear  system  are  the  vectors 


(  (x{yi)i/vf  \ 

/  a(xn+1)y/wn+1  \ 

<r(y2)V^ 

a{xn+2),Jwn+2 

Si  = 

and  S2  = 

\  <r(yi)Vvi  ) 

\  y/^n+m  / 

Note  that  the  scaled  values  of  a  charge  distribution  a  at  the  quadrature  nodes  yi,...,yi  are 
sufficient  to  fully  specify  the  restriction  of  the  <7  to  To-  Indeed,  a  coefficient  expansion  for  the 
restriction  in  terms  of  the  basis  (pi , ,(f>i  can  be  computed  by  applying  the  transformation  matrix 
constructed  in  the  preceding  section  to  the  vector  si. 

Remark  3.2.  The  procedure  of  this  section  applies  to  many  boundary  integral  operators,  not  just 
those  associated  with  Laplace’s  equation.  See,  for  instance,  [3]  for  a  discussion  of  the  application  of 
these  ideas  to  the  boundary  integral  operators  arising  from  the  Helmholtz  equation.  Moreover,  it  can 
extended  to  more  general  partial  differential  equations,  including  certain  classes  of  equations  with 
nonconstant  coefficients.  This  generalization  will  be  reported  at  a  later  date. 

Remark  3.3.  The  algorithm  of  [2]  produces  smaller  compressed  systems  than  the  procedure  of  this 
section,  but  it  does  so  at  the  cost  of  greater  complexity.  Rather  than  using  a  quadrature  for  the 
functions  (3.12)  to  discretize  the  interactions  of  the  curve  segment  To,  it  uses  a  quadrature  for  the 
smaller  collection  of  functions 
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The  advantage  of  this  approach  is,  of  course,  that  the  resulting  compressed  systems  of  linear  equations 
are  smaller.  The  disadvantage  is  that  interpolatory  quadratures  must  used  to  evaluate  the  integrals 

/  K(x,y)aj(y)ds(y) 

J  To 

when  x  £  9f2\To.  Depending  the  solver  applied  to  the  compressed  linear  systems,  this  approach  can 
actually  result  in  a  less  efficient  algorithm. 

4.  Numerical  Experiments 

This  section  describes  several  numerical  experiments  performed  to  assess  the  accuracy  and  stabil¬ 
ity  of  inner  product  preserving  Nystrom  discretizations  of  (1.4)  and  to  compare  such  discretizations 
to  those  obtained  via  the  standard  Nystrom  approach. 

Although  the  computations  described  here  could  be  accelerated  greatly  through  the  use  of  fast 
multipole  methods  and  fast  direct  solvers,  for  the  sake  of  simplicity  and  in  the  interests  of  repro¬ 
ducibility  the  experiments  were  conducted  using  standard  LAPACK  routines.  In  particular,  singular 
value  decompositions,  which  were  used  to  solve  all  linear  systems  save  for  those  related  to  the  com¬ 
pression  scheme  of  Section  3,  were  constructed  using  the  routine  DGESDD.  Linear  systems  arising 
from  the  compression  scheme  were  solved  using  the  DGESV  routine. 

Note  also  that  the  accuracy  of  solutions  was  given  priority  over  the  efficiency  of  discretizations. 
Substantially  more  efficient  discretizations  can  be  obtained  in  some  cases  by  sacrificing  a  few  digits 
of  precision. 

All  experiments  were  performed  on  a  PC  equipped  with  a  2.54GHz  Intel  Core  2  Duo  processor 
and  8GB  of  RAM.  Code  for  the  experiments  was  written  in  Fortran  77  and  compiled  with  the  Intel 
Fortran  Compiler,  version  11.1. 

4.1.  Graded  meshes.  Before  beginning  the  description  of  the  experiments  proper,  we  need  to 
introduce  the  terminology  which  will  be  used  throughout  this  section  to  describe  quadratures  for 
the  discretization  of  integral  operators.  Let  dCl  be  a  boundary  curve  parameterized  by  the  mapping 
r  :  [a,  b]  — >  dil  C  R2.  A  quadrature  formula  on  the  interval  [a,  b]  can  be  mapped  into  dil  via  the 
parameterization  r.  More  specifically,  the  quadrature  rule  X\, . . .  ,xn,W\  ,...,wn  on  [a,  b]  induces 
the  approximation  formula 


In  what  follows,  we  will  implicitly  identify  the  parameterization  domain  [a,  b]  and  the  boundary 
curve  dfl  through  the  use  of  this  mapping.  That  is,  we  will  describe  quadratures  on  the  curve  dfl 
by  describing  their  preimage  in  the  parameterization  domain  [a,  b]  under  this  mapping. 

To  that  end,  we  shall  call  the  piecewise  Legendre  quadrature  of  order  m  on  the  intervals 

[-2-j+1,-2~j],  [2-'\2-*+1],  j  =  l,...,N, 

a  simply-graded  mesh1  of  order  m  on  [— 1, 1]  with  cutoff  2~N .  We  call  the  image  of  such  a  quadrature 
under  the  substitution 

u  =  xk, 

where  k  is  a  positive  odd  integer,  a  graded  mesh  on  [—1,1]  of  order  to,  exponent  k,  and  cutoff  2~kN . 
Note  that  the  quadrature  formula  associated  with  a  graded  mesh  on  [—1,1]  is 


where  x\ , ...  ,xi  and  w\ ,...  ,wi  denote  the  nodes  and  weights  of  the  original  piecewise  Legendre 
quadrature. 

In  order  to  form  a  quadrature  for  a  boundary  curve  with  corner  points  parameterized  over  [a,  b) 
we  will  use  a  combination  of  graded  meshes  and  piecewise  Legendre  quadratures.  For  each  corner 

1The  term  simply-graded  is  borrowed  from  [7]. 


14 


JAMES  BREMER 


point  r{tj)  on  the  boundary  curve,  we  form  a  graded  mesh  on  [—1,1]  and  map  it  affinely  into  a 
neighborhood  [tj  —  S,tj  +  <5]  of  tj.  A  quadrature  on  the  remaining  portion  of  the  interval  [a,b\  is 
then  formed  using  a  piecewise  Legendre  rule.  In  many  of  the  following  experiments,  we  will  describe 
graded  meshes  using  these  4  parameters:  the  order  m  of  the  piecewise  Legendre  quadrature  used  to 
form  the  graded  mesh,  the  radius  S,  the  cutoff  ecut  for  the  mesh,  and  the  exponent  k. 

4.2.  A  family  of  domains  with  a  single  outward-pointing  corner.  For  0  <  6  <  7r,  let  f lg 

denote  the  domain  bounded  by  the  simply-connected  curve  parameterized  over  the  torus  [ — 7r,  7t]  by 

,  .  _  j—  2sin(f/2)  —7 r  <  t  <  0, 

6  |2sin(f/2)  0  <  t  <  tt 

yg(t)  =  —  tan(0/2)  sin(t). 

The  parameterization  (a ;g(t),yg(t))  has  counter-clockwise  orientation  and  maps  t  =  0  to  a  corner 
point  with  interior  angle  9  at  the  origin;  Figure  2  depicts  12^/2-  We  will  refer  to  the  fig  as  “snowcone” 
domains. 


Figure  2:  The  domain  fln/2- 


Figure  3:  The  single-layer  potentials  on  dfl^/2  (left)  and  (right)  which  generate  the  harmonic 

function  (4.1).  Note  that  the  function  on  the  right  has  an  unbounded  singularity  at  0  which  has 
been  truncated. 

Several  discretizations  of  the  operator 

Ta(x)  =  a(x)  +  [  ( log  \x  -  y\)  <r(y)ds(y) 

Jdn„/2  \ov x  ) 

were  formed  using  various  graded  mesh  quadratures.  In  each  case,  the  mesh  parameters  m  and  S 
were  30  and  1,  respectively,  and  an  180-point  piecewise  Legendre  quadrature  was  used  to  discretize 
the  operator  over  the  smooth  portions  of  the  curve  —  that  is,  the  regions  parameterized  over  the 
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intervals  [ — 7r,  — 1]  and  [1 ,  7r] .  For  each  graded  mesh  quadrature  considered,  both  the  inner  product 
preserving  discretization  (2.2)  and  the  usual  Nystrom  discretization  (2.3)  were  formed.  Moreover, 
each  inner  product  preserving  discretization  was  used  to  compute  a  solution  to  an  exterior  Neumann 
problem  on  the  domain  fln/ 2.  The  boundary  data  was  taken  to  be  the  restriction  to  dfln/ 2  of  the 
normal  derivative  of  the  function 

u(x)  =log|x-z0|  -log|a:-.2i|,  (4.1) 

where  zq  =  (1.0, 0.0)  and  z\  =  (1.1, 0.1).  Note  that 

u(x)  =  O  (l/|x|)  as  \x\  — >  00. 

For  each  obtained  solution,  the  error  in  the  representation  (1.3)  was  measured  at  300  points  on  the 
circle  C  of  radius  3  centered  at  the  origin.  The  results  appear  in  Table  1;  the  quantities  reported 
are  as  follows: 

•  ecut  refers  to  the  cutoff  for  the  graded  mesh; 

•  k  is  the  exponent  of  the  graded  mesh; 

•  n  is  the  dimension  of  the  systems  of  linear  equations; 

•  k  is  the  condition  number  of  the  inner  product  preserving  discretization; 

•  E  is  the  largest  absolute  error  observed  while  testing  the  representation  (1.3)  on  the  circle  C; 

•  njv  is  the  condition  number  of  the  standard  Nystrom  discretization. 
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Table  1:  Numerical  results  for  the  snowcone  domain  fln/2  ■ 


Figure  3  shows  the  single  layer  potential  on  dfln/2  which  gives  rise  to  the  harmonic  function  (4.1) 
in  f \c  .  It  is  clear  from  Table  1  that  inner  product  preserving  discretization  leads  to  well-conditioned 
systems  while  the  standard  Nystrom  approach  results  in  numerical  instability.  It  is  also  clear  that 
high  accuracy  can  be  achieved  for  this  problem. 

We  also  solved  exterior  Neumann  problems  on  a  collection  of  domains  O g  with  angles  9  close  to 
0  and  7 r  via  the  integral  equations 

a(x)  + ^  (^-log\x  -  y\j  a(y)ds(y)  =  2g(x),  x  e  dQ.e.  (4.2) 

For  these  experiments,  inner  product  preserving  discretizations  were  used  and  the  boundary  data  g 
was  taken  to  be  the  normal  derivative  of  the  potential 

log  \x  —  z2\  —  log  | a:  —  z3| ,  (4.3) 

where  z2  =  (1,0)  and  z 3  =  (1.1,0).  These  locations  were  chosen  because  z2  and  z 3  are  contained 
in  the  domains  fig  for  all  values  of  9.  Table  2  reports  the  results  of  these  experiments.  There,  9 
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refers  to  the  angle  of  the  corner,  n  is  the  number  of  quadrature  nodes  used  to  discretize  the  integral 
equation,  n  is  the  condition  number  of  the  system  of  linear  equations,  and  E  is  the  largest  absolute 
error  observed  while  testing  the  representation  (1.3)  at  a  collection  of  300  points  sampled  randomly 
in  the  box  with  corners  (—5,  —1)  and  (—4, 1). 
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Table  2:  Results  for  nearly  degenerate  snowcone  domains. 


It  is  expected  that  the  conditioning  of  the  linear  systems  discretizing  (4.2)  will  deteriorate  as  8 
goes  to  0  and  w.  This  phenomenon  is  not  related  to  the  corner  singularity.  It  is  instead  a  consequence 
of  the  fact  that  the  domains  Qg  become  elongated  as  8  goes  to  0  and  n.  The  same  behavior  can  be 
observed  in  potential  theoretic  operators  on  eccentric  ellipses.  Let  D  be  the  integral  operator 

DfM = * L  log  11  ■ vl) 

where  dE  is  the  counter-clockwise  oriented  boundary  of  the  ellipse 

©2+(t)2<1- 

The  operator  D  arises  in  the  solution  of  interior  Dirichlet  problems  for  Laplace’s  equation.  In 
particular,  I  +  D  plays  the  same  in  role  for  Dirichlet  problems  that  the  operator  I  +  N  does  for 
Neumann  problems.  We  can  assume  that  a  >  b  >  0.  Then  the  mapping  [—7 r,  7r]  — >  R2  defined  by 

x(t)  =  acos(f) 
y(t)  =  bsin(t) 


is  a  counter-clockwise  oriented  parameterization  for  dE.  It  allows  us  to  rewrite  D  as 

ab 


w-ii: 


-f(s)ds. 


(4.4) 


l_n  a2  +  b2  +  ( b 2  —  a2)  cos (s  +  t)  * 

The  spectrum  of  D  can  be  easily  computed  since  it  is  a  convolution  operator.  Its  eigenvalues 
{Anjwoo  are 

'a  — 


and 


An  — 


A— n,  — 


a  +  b 
a  —  b 


a  - 


n  =  0, . . . ,  oo, 


n  =  1, . . . ,  oo. 


It  follows  that  the  spectrum  of  the  operator  I  +  D  is  the  set 

'  a  —  b 


ct(7 +£>)  =  {  2}U  1± 


n  =  1, ....  oo 


From  this,  we  see  that  the  L2(dE)  condition  number  of  the  operator  I  +  D  is 


|/  +  D||2||(J  +  D)- 


=  1  + 


ON  THE  NYSTROM  DISCRETIZATION  OF  INTEGRAL  EQUATIONS  ON  PLANAR  CURVES  WITH  CORNERSI7 


This  condition  number  increases  as  the  ratio  of  a  to  b  becomes  large;  i.e. ,  as  the  ellipse  E  becomes 
more  eccentric.  A  similar,  but  slightly  more  complicated,  analysis  applies  to  the  operator  I  +  N 
associated  with  the  Neumann  problem  for  Laplace’s  equation. 

4.3.  A  family  of  domains  with  a  single  inward-pointing  corner.  In  this  section,  we  describe 
several  numerical  experiments  involving  a  family  of  “boomerang”  domains  which  have  a  single 
inward-pointing  corner.  For  ir  <  9  <  2ir ,  let  Tg  be  the  simply-connected  domain  whose  boundary 
is  parameterized  over  [71, 7r]  by 

.  .  f  2sin(3t/2)  —  7r  <  t  <  0, 

xg(t)  =  <  ,  ,  , 

^  —  2sin(3t/2)  0  <  t  <  ir 

y$(t)  =  tan(0/2)  sin(f). 

The  domain  T#  has  a  single  corner  point  with  interior  angle  9  at  the  origin;  see  Figure  4,  which 
depicts  \E'37r/2. 
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Table  3:  Numerical  results  for  the  boomerang  domain  'F37r//2- 


We  formed  several  discretizations  of  the  integral  equation 

a(x)+[  ( ^-log\x  -  y\)  a(y)ds(y)  =2g(x),  x€d^3v/2, 

Jdm 3„/2  \OVx  J 

with  g{x)  taken  to  be  the  normal  derivative  of  the  function  defined  by  (4.1).  A  number  of  different 
graded  meshes  were  used  and  both  inner  product  preserving  and  standard  Nystrom  discretizations 
were  formed  for  each  mesh.  In  all  cases,  the  parameters  S  and  m  were  taken  to  be  1/2  and  30, 
respectively,  and  the  integral  equation  was  discretized  over  smooth  portions  of  the  curve  using 
piecewise  Legendre  quadratures.  The  results  appear  in  Table  3;  the  quantities  reported  are  the 
same  as  those  in  Table  1  of  the  preceding  section.  The  charge  distribution  a  which  gives  rise  to  the 
harmonic  function  (4.1)  in  tUg  ,2  is  shown  in  Figure  3.  Note  that  a  has  an  unbounded  singularity 
at  0. 
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Figure  5:  The  domain  'F27r_7r/io- 


We  also  solved  exterior  Neumann  problems  on  a  collection  of  domains  Tg  with  angles  9  close  to 
7 t  and  27t  via  inner  product  preserving  discretizations  of  the  integral  equations 

<x{x)  +  j  (^J^\og\x-y\SJa(y)ds(y)  =  2g(x),  x  £  d'&g.  (4.5) 

For  each  problem,  the  boundary  data  was  taken  to  be  the  normal  derivative  of  the  function  (4.3). 
Table  4  presents  the  results  of  these  experiments;  the  quantities  reported  are  the  same  as  those  in 
Table  2.  Note  that  the  resulting  charge  distributions  exhibit  unbounded  singularities  in  the  case  of 
angles  near  i r  but  not  in  the  case  of  angles  near  27t;  see,  for  instance,  Figure  6,  which  shows  the 
single-layer  potentials  on  94'2w_w/io  and  S\E',r+7r/10  which  generate  the  harmonic  function  (4.3). 


Figure  6:  The  single-layer  potentials  on  94'2w_7r/10  (left)  and  d'^n+n/w  (right)  which  generate  the 
function  (4.3). 


Nonetheless,  much  more  difficulty  was  encountered  in  the  solution  of  the  Neumann  problems 
when  9  was  close  to  2n  than  when  9  was  close  to  7r.  These  difficulties  do  not  arise  because  of  the 
corner  singularities,  but  because  of  the  elongated  segments  of  the  boundary  contours.  Figure  5 
shows  the  domain  4'27r_W//10,  whose  boundary  features  two  elongated  regions. 


9 

n 

K 

E 

7 r  +  7r/10 

1680 

7.59x10+° 

l.llxlO"14 

7T  +  7t/20 

1800 

1.22x  10+1 

7.55xl0"15 

7T  +  7T/100 

1860 

5.53xl0+1 

5.40xl0"14 

7T  +  7T/1000 

4080 

5.56xl0+2 

3.70xl0-12 

27T  —  7r/10 

4440 

2.09xl0+3 

2.46xl0"14 

to 

1 

to 

o 

2280 

1.06xl0+4 

4.33xl0-14 

2n  —  7r/100 

3808 

2.85xl0+5 

6.84xl0-14 

Table  4:  Results  for  nearly  degenerate  boomerang  domains. 


4.4.  A  domain  with  8  corner  points.  In  this  experiment,  we  solved  an  exterior  Neumann  problem 
on  the  simply-connected  “inkblot”  domain  Oink  shown  in  Figure  7.  The  boundary  of  this  domain 
has  8  corner  points,  4  inward-pointing  and  4  outward-pointing,  and  is  parameterized  by  the  polar 
equation 


r(9)  =  4  +  2  |cos  (40) |  sin  (4 9) ,  0  <  9  <  2n. 
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The  boundary  data  g  was  taken  to  be  the  normal  derivative  of  a  potential  generated  by  two  charges 
of  equal  magnitude  and  opposite  signs  which  were  randomly  placed  in  the  interior  of  flink. 


Figure  7:  The  domain  Qink  of  Section  4.4. 


n 

^corner 

^"smooth 

K 

E 

?SVD 

T 

-L  compress 

Uncompressed  system 

8700 

1020 

540 

1.95xl0+1 

9.27xl0"15 

1112.5 

- 

Compressed  system 

1004 

58 

540 

5.41xl0+1 

5.55xl0_ls 

1.34 

3.26 

Table  5:  Numerical  results  for  the  domain  flink  of  Section  4.4.  Entries  marked  with  a  dash  are  not 
applicable. 


In  order  to  obtain  the  solution  of  this  Neumann  problem,  an  inner  product  preserving  discretiza¬ 
tion  of  the  integral  equation 

cr(x)+  f  ( ^-\og\x  -  y\)  a(y)ds(y)  =  2g(x),  x  €  dflink, 

Jdnink  \(Jvx  ) 

was  formed  and  a  singular  value  decomposition  of  the  matrix  associated  with  this  linear  system  was 
computed  and  used  to  solve  the  linear  system.  The  resulting  charge  distribution  is  singular  at  each 
corner  point  of  dflink;  the  singularities  occurring  at  inward-pointing  corner  points  are  unbounded 
while  those  at  outward-pointing  corner  points  are  bounded.  The  compression  scheme  of  Section  3 
was  also  applied  to  this  discretization  in  order  to  form  a  compressed  system  of  equations.  Table  5 
presents  the  results;  the  quantities  reported  are  as  follows: 

•  n  is  the  dimension  of  the  system  of  linear  equations; 

•  ficomer  is  the  average  number  of  nodes  used  to  discretize  corner  regions; 

•  nSmooth  is  the  number  of  quadrature  nodes  used  to  discretize  the  smooth  portions  of  the  curve; 

•  k  is  the  condition  number  of  system; 

•  E  is  the  largest  absolute  error  in  the  representation  (1.3)  occurring  on  the  circle  of  radius  6  centered 
at  the  origin; 

•  Tsv d  is  the  wall-clock  time,  in  seconds,  which  was  required  to  compute  the  singular  value  decom¬ 
position  of  the  matrix; 

•  TcompTess  is  the  wall-clock  time,  in  seconds,  required  to  execute  the  corner  compression  procedure. 
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4.5.  A  polygonal  domain  with  38  corner  points.  In  this  final  experiment,  an  inner  product 
preserving  discretization  of  the  the  integral  equation 

a(x)+[  ( -^-log\x-y\ )  a(y)ds(y)  =  2g(x),  x  €  <9Ustar, 

•'3Qstar  \°vx  J 

where  Qstar  is  the  simply-connected  “starburst”  domain  shown  in  Figure  8,  was  formed.  The  bound¬ 
ary  of  fl3tar  is  a  polygon  with  38  vertices.  The  procedure  of  Section  3,  which  is  absolutely  essential 
in  this  case,  was  applied  to  that  discretization  in  order  to  form  a  compressed  system  of  linear  equa¬ 
tions.  The  compressed  linear  system  was  used  to  solve  the  exterior  Neumann  problem  on  f2star  for 
boundary  data  which  was  the  normal  derivative  of  a  potential  generated  by  two  charges  of  equal 
magnitude  and  opposite  sign  placed  randomly  in  the  domain.  Table  6  describes  the  results;  the 
quantities  reported  are  the  same  as  those  of  Table  5  in  Section  4.4. 


Figure  8:  The  starburst  domain  Qstar 


n 

^corner 

^smooth 

K 

E 

Tsvo 

T 

-L  compress 

Uncompressed  system 

70740 

1800 

2340 

* 

* 

* 

- 

Compressed  system 

5198 

75.21 

2340 

1.40xl0+2 

1.45x  10-15 

223.2 

84.2 

Table  6:  Numerical  results  for  the  domain  fl3tar.  Entries  marked  with  a  dash  are  not  applicable, 
while  those  marked  with  an  asterisk  could  not  be  calculated  due  to  computational  limitations. 


5.  Conclusions 

There  is  an  extant  problem  regarding  the  solution  of  boundary  integral  equations  on  planar 
domains  with  corners,  but  it  is  not  the  stability  or  accuracy  of  graded  mesh  discretization.  Rather, 
the  outstanding  issue  is  the  compression  of  the  large  systems  of  linear  equations  which  result  from  the 
discretization  of  boundary  integral  equations  on  such  domains.  Of  particular  interest  are  schemes, 
like  those  presented  in  [7]  and  [3,  2],  which  do  not  require  a  priori  analytic  information  about  the 
behavior  of  solutions.  This  property  allows  them  to  be  applied  to  the  many  problems  which  can 
be  formulated  as  boundary  integral  equations  with  little  modification.  Moreover,  it  is  almost  a 
prerequisite  for  schemes  in  the  three-dimensional  environment,  where  estimates  of  the  behavior  of 
solutions  are  difficult  to  obtain.  Indeed,  one  of  the  principal  advantages  of  the  approach  of  Section  3 
is  that  it  admits  generalization  to  three-dimensional  domains  with  singularities,  a  generalization 
which  will  be  reported  by  the  author  at  a  later  date. 
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